Non-local exchange effects in zigzag edge magnetism of neutral graphene nanoribbons 
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We study the role of non-locality of exchange in a neutral zigzag graphene nanoribbon within the s-orbital 
unrestricted Hartree-Fock approximation. Within this theory we find that the magnetic features are further 
stabilized for both the intra-edge and inter-edge exchange compared to mean field theories of zigzag ribbon 
based on local exchange (like the Hubbard model or the ab-initio Local Density Approximation). The inter- 
edge exchange produces an enhancement of the band-gap of the magnetic ground-state solutions. The effect 
of this enhanced exchange on the edge states can not be satisfactorily achieved by a local interaction with 
renormalized parameters. 

PACS numbers: 75.30.Et, 75.75.+a, 73.22.-f, 73.20.-r 



I. INTRODUCTION 



Zigzag terminated graphene ribbons have recently at- 
tracted a lot of attention both by theorists^ " 23 ' 25 and 
experimentalists 26 " 39 partly due to the surge in popularity 
of graphene after seminal transport experiments^!^ 2 - The 
unique properties of zigzag edge localized states has prompted 
using zigzag terminated nanoribbons as a testbed for propos- 
als of magnetism in graphitic systems 43 or for illustrating ex- 
otic physics related to the bulk topology of single sheet of 
multilayer graphene in presence of strong spin-orbit coupling 
or magnetic fields 43rd 7 . Recent experiments have found traces 
of localized spin polarization at zigzag edges of graphen o 37 " 39 
in agreement with theories predicting spin polarization at 
zigzag edges in presence of electron interactions ! 1 i 7 i 9 i 11 The 
agreement of the qualitative features between solutions pre- 
dicted by different methods such as the Hubbard model and 
DFT with different semi-local approximations for the en- 
ergy functionals point towards a robust underlying property of 
these systems for which the magnetic solutions are not very 
sensitive to the specific details of electron interaction. Con- 
sidering that magnetic ordering is a long-ranged phenomena 
one natural question that follows is how much physics we are 
missing by an inadequate treatment of non-local exchange in 
our solutions. The Hubbard model used in several previous 
papers 13 ' 5Q i 54 ' 55 relies on a strictly local onsite interaction term 
while the common approximations of the energy functionals 
in DF T 56 ' 57 are modeled from local or semi-local exchange 
and correlation energies 58 of the homogeneous electron gas. 
In both cases the description of the exchange hole tends to 
be excessively short ranged in addition to uncertainties asso- 
ciated with self-interaction errors in the case of DFT. In this 
work we revisit the problem of edge magnetization, predicted 
by mean field theories in zigzag nanoribbons, within Unre- 
stricted Hartree Fock (UHF) approximation with non-local 
exchange, that does not suffer from the shortcomings of the 
previous approximations, and study systematically different 
truncation ranges for the electron interaction. 



We find that a non-local interaction enhances the gaps in 
the band structure, a feature that can be traced back to both 
the inter-edge tunneling term and intra-edge magnetic order. 
As we discuss below models with short-range interaction are 
unable to mimic the large size of the gap in a consistent way 



although the solutions remain similar in many aspects. 

We start in section II introducing the Hartree-Fock theory 
in zigzag ribbons and cast it in the form of a two dimensional 
edge state bands mode l 13 ' 15 extending beyond the formulation 
based on the Hubbard hamiltonian. In section III we show 
within this framework how the non-local exchange modifies 
the effective hamiltonian matrix elements when we extend the 
onsite interaction Hubbard model to include farther neighbor 
interaction terms. We then move on to discuss in section IV 
the localization properties of the edge state wave functions 
along the direction of the ribbon discussing the dependence of 
the intra-edge band-gap as a function of the electron interac- 
tion strength. Subsequently in section V we assess the impact 
of non-local exchange on the energetics of the inter-edge an- 
tiferromagnetic ground state. We finally close the paper with 
a summary and conclusions section. 



n. TWO BANDS HARTREE-FOCK THEORY IN ZIGZAG 
RIBBONS 



A simple and yet fairly accurate way to study edge mag- 
netism in a zig-zag ribbon in a mean field theory is to calcu- 
late, for each Appoint, the effect of interaction on the edge state 
wave functions not allowing edge and bulk to mix their wave 
function s 13 ' 15 . These edge state wave functions (see appendix 
for more details) become exponentially localized at the edges 
whenever 2n/3a + 1 /W < \k\ for sufficiently wide ribbons 5 - 60 
where the ribbon width is W = V5aN/2 and N is the number 
of atom pairs in the unit cell across the ribbon and the lattice 
constant is a = 2.46A. We label in an abbreviated notation the 
ribbons edge states as \k—) and \k+), antisymmetric and sym- 
metric wave functions across the ribbon, that can be symmet- 
rically and antisymmetrically combined to obtain basis func- 
tions that are mostly centered either at the (L)eft or (R)ight 
edge in the ribbonji 5 - We denote the respective ^-dependent 
wave function amplitude at each lattice site I in the unit cell 
as and R^. We can use this new basis to represent the 
Hamiltonians as a two by two matrix for each spin rj =t / I 



H„(k) 



H L L,o(k) Hlr, a (k) 
Hrl, a (k) H RR:(7 (k) 



(1) 
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The tight-binding band term Hamiltonian is described by a 
simple nearest neighbor hopping term 70 = — 2.6eV. In this 
LR basis representation the Hamiltonian would consist of a 
tunneling term mixing states located at both edges in the rib- 
bon for each spin 



H TB (k) = 



t TB (k) 
t TB (k) 



(2) 



where tjB (k) is the energy dispersion for the tight-binding 
edge conduction band ( i.e., the basis we use consists of two 
eigenstates of H TB (k)). For a neutral ribbon the electrostatic 
Hartree term cancels with the positive background charge and 
contributes only with a constant term Co in the diagonal ele- 
ments that can be chosen at our convenience. Then the inter- 
action term of the Hartree-Fock hamitonian projected on the 
two-bands LR basis effectively reduces to an exchange contri- 
bution on both the diagonal and off diagonal matrix elements 



H LL . a {k) 

H LR .o{k) 
H RR , a {k) 
where we have defined 



Co - ^L k iLki'Fx\ a (k) 

W 

t TB (k)-Y,L kl R kl >F^ a (k) 
IV 



Fx,c{k) = L U x {k'-k)(n k ,„,. a ) 

k! 



(3) 
(4) 
(5) 

(6) 



where 



n k!W. 



is the density matrix in the lattice Bloch 



states and U% (k' — k) represents exchange Coulomb integrals 
whose explicit expressions are presented in the appendix. The 
Coulomb integrals between the 7T-orbitals located at different 
sites can be approximated in Hartree atomic units by effective 
terms of the form 64 -^ 



Vai{d) 



1 



£ r y/al + d 2 



(7) 



where d is the distance between lattice sites and the bond- 
ing radius of the carbon atoms a = a/ (2v3) account for a 
small damping of the interaction due to the shape and finite 
spreading of the wave function around the lattice. The on- 
site repulsion term remains unknown and is commonly brack- 
eted between U = 2eV and U = 6eV4^ In our calculation 
we used a relatively small value of U = 2.5eV, considering 
that the onsite repulsion strength of U = 2eV in the Hubbard 
model where the Coulomb tail is neglected is enough to re- 
produce the LDA band structures 13 . We include in our calcu- 
lations an effective screening of e r — 4 resulting in a global 
damping of the interaction strength by a factor 1/4 to account 
for possible effects of dielectric screening (usually e r = 2.5 
in S1O2) and additional screening effects due to sp% orbitals 
that are neglected in our approximation. The range of the 
Coulomb interaction is truncated to up 16 nearest neighbors. 
With the above prescriptions we obtain for the ground state of 
the neutral system an intra-edge interaction gap near ka ~ n 
of 1 AeV midway between the B3LYP prediction^ of 2.2eV 



and the LDA prediction 7 of 0.5eV. We write the two-bands 
effective Hamiltonian in the following way: 

H a (k) = A AF (k)ax z +A F (k)aI+t a (k)x x , (8) 

where the pauli matrices Xn and / are acting in the space of left 
and right edges and a is +(— ) for up(down) spin. The renor- 
malized inter-edge tunneling t a (k) is given by the bare hop- 
ping plus a spin dependent enhancement from the exchange 
interaction. The intra-edge exchange potentials A AF (k) and 
A F (k) correspond to a particular choice of A AF ^ (k) or A F ^ (k) 
such that they are positive quantities at k ~ nja. Written in 
terms of the matrix elements defined in Eq. ( 131415b for each 
momentum k they are: 



ta (k) 
Af(k) 



HlR, a (k) 

1 



2 (H L L,a{k)-H RR ^{k)) 
&o(k) = \{H L L,a+H RR . a ){k). 



(9) 
(10) 

(11) 



Note that we have assumed a collinear spin arrangement and 
there is no spin-mixing term. Therefore the 4x4 Hamiltonian 
matrix is reduced to two block diagonal sub-matrices in Eq. 
((H) for each spin sub-space. The energy dispersions associ- 
ated with the edge states of the two-band Hamiltonian can be 
written as 



E±(k) = oA*(k)±^A AF '(k)+tl(k) 



(12) 



In the AF case the term A F (k) in equation ([8]) amounts to an 
edge-independent value for both up and down spin and can 
be set to zero. We notice that the exchange spin splitting 
terms given by the diagonal elements A AF (k) have opposite 
signs for left and right edge states, and a global sign rever- 
sal when we consider the opposite spin Hamiltonian. On the 
other hand, in ferromagnetic self-consistent solutions we get a 
nonzero net spin polarization in the sample due to edge elec- 
trons with spins pointing in the same direction. In this case be- 
cause A AF = the intra-edge exchange energy gain included 
in A F (k) is reflected in the overall shift of the band energy 
for each k. The opposite spin state is shifted in an equal and 
opposite direction and the bands split in an amount of 2A F (k) 
at each £-point. 



III. NON-LOCAL EXCHANGE AND THE EFFECTIVE 
HAMILTONIAN MATRIX ELEMENTS 

The non-locality in exchange appear when the Coulomb in- 
teraction range is extended beyond the onsite repulsion and 
is most easily explored comparing the Hartree-Fock Hamil- 
tonian matrix elements of equations ( 131415b with those we 
obtain using the Hubbard onsite repulsion model. We have 
used the value of U = 2eV for our reference Hubbard model 
calculations following the convention of previous work a 13 ' 15 
to match the LDA band gaps. In Fig. 1 we show a com- 
parison of the band structures in the HF and Hubbard model 
for the lower energy AF and higher energy F configurations. 
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FIG. 1 : (Color online) Comparison of HF and Hubbard model band 
structures for the lower energy AF (left panel) and meta-stable higher 
energy F (right panel) spin configurations for a nanoribbon with N = 
20 carbon atom pairs in the unit cell. In the HF solution we find 
enhanced band-gap openings due to the increase in strength of both 
intra-edge exchange manifested near ka^n and inter-edge tunneling 
that become more significant as we move away from this point. 



The effects of non-local exchange in the band structure can 
be understood more clearly separating their contributions for 
the intra-edge exchange potentials corresponding to the diag- 
onal elements of the effective Hamiltonian obtained through 
equations ( 13151101 ) and the off diagonal inter-edge tunneling 
given in equations d4|91 l. We discuss our comparison only for 
AF matrix elements because the effective matrix elements for 
F solutions remain qualitatively identical. These matrix ele- 
ments are represented in Fig. (2) for different choices in the 
truncation of the interaction ranges for the Coulomb integrals 
in equation @. Since the left and right basis functions are 
strictly localized on one of the sublattices the contributions 
to the diagonal term A (k) are due to interaction couplings 
within the same sublattice, therefore these terms contain in- 
formation of exchange within the same edge atoms. One in- 
teresting feature we observe is that the function A AF (k) ob- 
tained for the fully non-local calculation can be effectively 
reproduced in most of the edge states zone of the Brillouin 
zone (\k\ < 2n/3a + 1 /W) by the Hubbard Hamiltonian with 
a larger effective onsite repulsion. This effective repulsion 
takes the value of f/ e ff = 5.3 eV for the interaction param- 
eters we have chosen for the HF calculation, which is very 
close to the critical value of £/c/l7o| = 2.23 above which 
the choice would be unphysical because the ground-state of 
a 2D graphene sheet develops an antiferromagnetic spin den- 
sity wave solution!. For every choice in the cutoff for the 
Coulomb interaction range we could obtain different effective 
choices of t/ e ff that can reproduce A AF (k) satisfactorily. In 
the Hubbard model the k dependence of the diagonal terms of 
the Hamiltonian in Eqs. (|3] 01 is manifestly due to the coef- 
ficients Li, of the basis function with the largest contributions 
coming from the lattice sites at the edge and decaying expo- 
nentially as we move into the bulk. In the corrections due to 
farther neighbor interactions (i.e., beyond the Hubbard model) 
still the dominant contribution to the matrix element A for 
a given distance of electron interaction are those connecting 
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FIG. 2: (Color online) Enhancement of the matrix elements of the 
effective hamiltonian defined in equations d9l 1 01 lib evaluated in the 
HF approximation on a ribbon with N = 20. The subscripts in the la- 
bel indicate the number of nearest neighbors considered in the trun- 
cated range of the effective Coulomb interaction. We do not repre- 
sent the results for F configuration because the behavior is similar to 
what is found for AF solutions. Left Panel: Values of the intra-edge 
exchange potential A AF (k) obtained in the HF approximation for dif- 
ferent interaction ranges. With an appropriate choice of an effective 
U e ff = 5.3eV in the Hubbard model we can match A AF (tz/o) of the 
HF calculation for most of the region in the Brillouin zone where 
the wave functions are edge localized. The red curve represents the 
effective Hubbard model that agrees well with the HF calculation in 
almost the entire range of edge states zone \k\ < 2n/3a+l /W. Right 
Panel: Enhancement of tunneling amplitudes tf n {k) for AF solu- 
tions in presence of long-ranged interaction representing the greater 
coherence between wave functions located at each sublattice, and 
therefore each edge in the ribbon. For onsite interactions the tunnel- 
ing has the same form as the tight-binding conduction band disper- 
sion. 



to the edge ribbon atom, resulting in a ^-dependent behavior 
proportional to that of the h\ located at the edge sites. This 
remarkable behavior for which every additional neighbor con- 
tribution in the exchange potential has the same fc-dependent 
behavior in the edge states region is most likely the reason 
for the good overall agreement of the edge state properties in 
zigzag ribbons calculated within a simple Hubbard model and 
other mean field calculations. Outside this region of the Bril- 
louin zone we find substantial differences of the potentials for 
the Hubbard model and Hartree-Fock hamiltonian matrix ele- 
ments but these differences do not introduce relevant changes 
in the magnetic configuration of the edges. From the results 
shown above we expect that the quantitative agreement of the 
gap at ka ~ n between LDA approximation and the Hubbard 
calculations 13 with U = 2eV suggests that the LDA will have 
a tendency to underestimate the energetics of the ferromag- 
netic spin alignment along a zig-zag edge. In fact the LDA 
approximation 56 would introduce a larger gap opening if the 
spin density of the core electrons of carbon were considered 
when evaluating the spin dependent LDA exchange potential. 

The off diagonal tunneling term t a (k) given in equation 
(0 consists of a tight-binding term and an interaction term 
coupling different sublattices. The onsite interaction term in 
the Hubbard model cannot couple different sublattices, and 
therefore the two edge sites, and the matrix elements reduces 
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to the tight-binding band dispersion of the conduction edge 
band 15 . This matrix element can be enhanced in the presence 
of the long-ranged interaction. The enhanced coherence be- 
tween left and right edge solutions lead to an increase in the 
band-gap between the valence and conduction bands of the 
antiferromagnetic solution. This enhancement of the tunnel- 
ing term t a (k) is more pronounced for ^-points outside the 
edge states zone because the density tails spread more into the 
bulk. However it remains essentially zero in the surroundings 
of ka ~ % since the wave functions are strongly localized at 
the ribbon edge borders. 

The ribbon width dependence of both Hartree-Fock and 
Hubbard solutions remain similar. Previously 15 we had illus- 
trated using the Hubbard model that the ribbon width depen- 
dence of the solutions can be summarized through the behav- 
ior of the Hamiltonian matrix elements near the valley point 
k ~ 2n/3a that obey the following scaling rules as a function 
of ribbon width W 

A AF l F {k) = W- l A AF/F (qW) (13) 
t(k) = §t(qW) (14) 

where q = ka — 2n/3 is the dimensionless momentum mea- 
sured from the valley point and the functions A AF I F , t are 
functions that do not depend on the width of the ribbon. These 
rules for wide enough ribbons are not altered with the pres- 
ence of non-local exchange shown in equations ([3] |U [5]) be- 
cause the additional non-local terms are quadratic in and 
Rki- These in turn follow a similar scaling rule that can be 
obtained using a continuum model^ 

Rki = W-^ 2 R qW>l (15) 
L kl = W-^ 2 L qWJ (16) 

where R q w i and L q w / are scaled functions that do not depend 
on the width of the ribbon. Hence the band-gap size and the 
position of the gap as well as the up and down spin crossover 
point of the F solution (measured from the valley point) all 
follow a W~ ribbon width dependence. 



IV. INTRA-EDGE LOCALIZATION PROPERTIES OF 
EDGE STATES 

We explore here the properties of localization of the zigzag 
edge states along the x axis, the direction of periodicity in the 
ribbon. We estimate this quantity in real space evaluating a 
partial Fourier transfom of the £-point dependent density com- 
ponent centered at lattice site I in the ribbon unit cell 

F la {x)= ( dke- ikx \V kl {v)\ 2 (17) 

Jk>\2x/3a+l/W\ 

where (r) is the I th component of edge band Bloch wave 
function. The integration in k space is over the region of edge 
states which give the relevant contribution to the edge spin 
polarization. This definition, reminiscent of the way Wannier 
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FIG. 3: (Color online) Wave function localization along the rib- 
bon centered at different sites in the ribbon unit cell calculated for a 
ribbon with N = 12 corresponding to the antiferromagnetic ground 
state solution and interaction free tight-binding solution. We can no- 
tice a substantial spreading of the wave function at the ribbon edge 
enveloping about 7 units of lattice constant a along the ribbon. 



functions are calculated, shows us the way the electrons lo- 
calize in real space along the edge. In figure (f3j) we focus on 
the row of atoms at the ribbon edge and we show that the edge 
state wave functions spread about seven atomic lattice sites 
on average, and this spreading length is essentially same for 
paramagnetic tight-binding solutions and the solutions with 
edge spin polarization. These features are robust to changes 
in the ribbon width and the details of the Coulomb interac- 
tion with negligibly small departures from this form when 
we change the strength of the Coulomb interaction or modify 
the range of the Coulomb interaction. This robust feature of 
the wide smearing of the electron density on the edge atoms 
into the neighboring atoms located several lattice constants 
away along the edge direction is key in producing the long 
reaching ferromagnetic spin correlation length in a zigzag rib- 
bon edgeM^ Due to this far reaching overlap of the electron 
density along the edge even short ranged exchange interac- 
tions are able to connect edge states centered around different 
edge atoms and induce an energetically favorable parallel spin 
alignment. 

An intuitive way to relate this edge wave function local- 
ization length along the edge with the band structure or the 
ribbons near ka ~ n is by relating the gain of exchange en- 
ergy per particle for the occupied electrons near the edge 
atoms through A max — A AF I F {it/a) with an effective local- 
ization length A through the formula 

e 2 

&max — 7T- (18) 

e r A 

We present in table Q] the values of A max and X for differ- 
ent choices of £, calculated for a ribbon of N = 12 atom 
pairs width. In this analysis we define the onsite term as 
U =\0/e r eV making it dependent of the dielectric screening. 
The value of A max remain practically constant as we make the 
ribbon wider and already for N = 12 it gives a good estimate 
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of the infinite width limit. The resulting values of X shows 





1 2 3 


4 


5 


^max 


2.35 1.16 0.78 


0.58 


0.46 


X 


2.49 2.52 2.5 


2.5 


2.5 



TABLE I: Effective localization length X (in units of lattice con- 
stant a = 2.46A) estimated from the shift from tight-binding bands 
k-max (in eV) near ka ~ n for antiferromagnetic spin configuration 
in graphene nanoribbons. We can observe that X does not change 
much as a function of the values of dielectric constants considered. 
We can observe that the magnitudes of X are consistent with the es- 
timations that can be extracted for the localization lengths presented 
in the previous section. 

only a very small variation when different values of the rela- 
tive dielectric constant e r are used, in agreement with the fact 
that the shapes of the density localization remains practically 
unchanged. 



V. ENERGETICS OF INTER-EDGE EXCHANGE 
COUPLING 

The nature of antiferromagnetic spin polarization of edge 
states on different sublattices (and therefore different edges) 
was studied previously as an unusual type of superexchange 
interaction 15 where it is energetically favored with respect to 
the ferromagnetic state both in the kinetic energy and interac- 
tion energy. For the Hubbard model the inter-edge coupling 
is mediated by the band Hamiltonian tunneling and the local 
spin-polarization is provided by intra-edge exchange. In the 
presence of non-local interactions the inter-edge off-diagonal 
tunneling term plays a role in determining the total exchange 
energy. In this section we examine the impact of non-locality 
of exchange in the energetics of the ground-state solution by 
obtaining the differences in the total energies between AF and 
F mean field solutions. In a neutral nanoribbon the electro- 
static Hartree energy is the same in both AF and F configura- 
tions and the total energy difference per edge carbon atom AE 
consists of the kinetic energy term and the exchange term 



AE = E F - E AF = AT + AE X . 



(19) 



The kinetic energy difference can be calculated from one- 
body averages of the tight-binding Hamiltonian for each oc- 
cupied edge-band states and is essentially the same derivation 
as in reference B15I1 while for the exchange energy we need to 
sum two body exchange integrals. It will be useful to write 
it as an integral of a ^-dependent function defining implicitly 
£x(k) 



1 occ f 
E x = -iL K U= / dke x (k) 

J BZ 



(20) 



where Kjj is the standard definition of exchange integral as 
can be found in reference 161]. The labels i, j represent occu- 
pied single particle states including the k quantum numbers, 



band labels and spin index. In Fig. (4) we show the exchange 
energy difference as a function of Appoint in the Brillouin zone 
for the Hubbard model and the Hartree-Fock approximation. 
We can observe that the details of fc-point dependent contribu- 
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FIG. 4: (Color online) Comparison of fc-point resolved exchange en- 
ergy differences between Hubbard and HF calculations. Left panel: 
Exchange energy difference density Ae^ (k) for Hubbard calculation 
with U = 2eV Most of the exchange energy difference comes from 
the regions of \k\ ~ 27r/3a while the region ka ~ n also adds a small 
contribution. Right panel: In presence of longer ranged interactions 
we find contributions to the energy differences in a wider range of 
fc-points around the valley point 2k /3a thanks to the enhanced intra- 
edge tunneling terms. 



tions to the exchange energy difference is substantially modi- 
fied in presence of non-local interactions. The global enhance- 
ment in magnitude leads to a further stability of AF solutions. 
In particular the off diagonal interaction terms in the effec- 
tive Hamiltonian further enhances the stability of AF solutions 
through an enhanced inter-edge tunneling. 

The ribbon width dependence of the total energy differ- 
ences can be derived from similar considerations as in ref- 
erence [15] and follows a W~ 2 decay law. 



VI. SUMMARY AND CONCLUSIONS 

In this paper we have addressed the effects of non-locality 
of exchange in the edge states of a neutral zig-zag graphene 
ribbon from different points of view. We started with a for- 
mal introduction on how the far reaching interaction terms 
between distant sites can influence the values of the matrix 
elements of the effective two-bands Hamiltonian describing 
the edge states, distinguishing the different roles in terms of 
intra-edge and inter-edge interaction manifested respectively 
in the diagonal and off diagonal matrix elements, representing 
respectively the direct exchange energy giving rise to the spin 
polarization and the tunneling between both sublattices that 
mix states localized at different edges. Further insight was 
gained studying the properties of localization of edge states 
along the direction of the ribbon finding for those states sitting 
on the edge atoms a rather long ranged enveloping function in 
the direction of periodicity, a fact that would make possible 
even for very short ranged interaction terms to connect with 
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FIG. 5: (Color online) Total energy difference per edge carbon atom 
as a function of ribbon width showing a W~ 2 decay behavior. The 
fitting curves in units of eV versus A are 6/(W 2 + 300) for AE^f , 
1 /(W 2 + 50) for AE^f with U = 2eV and 2.6/ (W 2 - 100) for U = 
53eV. 



to those obtained by more local prescriptions like LDA/GGA 
functionalsii 
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the wave functions centered several lattice constants away in 
the direction of the ribbon. Finally we analyzed the impact of 
non-local interaction in the energetics of inter-edge coupling 
largely responsible for the antiferromagnetic spin polarization 
of the system. 

In light of our studies we have been able to understand bet- 
ter why the Hubbard model based on a strictly short ranged 
onsite interaction has turned out to give results consistent to 
solutions found in more elaborate studies. The main reasons 
would be on the one hand that for edge states the inter-edge 
tunneling not captured by the Hubbard model is relatively less 
important with respect to the intra-edge spin splitting contri- 
butions whose functional form in £-space is accurately cap- 
tured by an onsite interaction. On the other hand, the spread- 
ing of the edge state wave functions along the ribbon allows 
edge states centered at different atomic sites to connect each 
other even for strictly short ranged onsite interactions. Edge 
width depending scaling of band gaps and energy differences 
of different spin polarized states do also follow a similar law 
of and W~ 2 respectively as found in the Hubbard model 
since the width dependence of those quantities are dictated by 
the wide ribbon asymptotic behavior of the edge state wave 
functions. Even though the above mentioned qualitative fea- 
tures remain the same both for non-local and local exchange 
described by the HF and Hubbard models we have also shown 
that the discrepancies in the total energy differences and the 
exchange potentials outside the edge states region of the Bril- 
louin zone cannot be modeled accurately with a renormalized 
effective onsite interaction term. This effect can have a rel- 
evant influence in calculating the the spin stiffnes a 12 ' 16 or in 
the band structure of the system when it is shifted away from 
the neutrality poin t 20,21 The analysis we carried out illustrates 
the role of non-local exchange in altering the band dispersion 
and band-gap size due to spin polarized edge states, explain- 
ing the relative enhancement of band gaps in the results ob- 
tained with non-local B3LYP 9 type functionals with respect 
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Appendix I. Hartree-Fock formulation in the two band analysis where \]/ kX (r) are the Bloch state wave functions. 



The full Hartree-Fock Hamiltonian can be reduced to a two 
dimensional matrix representation for each spin. We will rep- 
resent this Hamiltonian matrix in a suitably chosen basis func- 
tion, for example the LR basis we have defined in the main 
text. Given the basis of Bloch functions 

(r|U) = y/u (r) = JL£ (r - R, - z,) rj CT 



(21) 

where 77^ represents the spinor, T/ represents the displacement 
vector for sublattices in the unit cell that can be labeled with 
/, and Nk is the total number of ^-points, or equivalently, the 
number of unit cells in the system repeated in the periodic 
direction. The label A = (I, a) represents both the lattice site 
label I and the spin a. In this basis the general expression of 
the Hamiltonian is 



where 



and 



kXl' 



c kX c kX 



£ U% x ' (k! - k) U, x ,cy X ) c\ x c kX , (22) 

k'XX' 



(kAk'A'|V|kAk'A') (23) 
dndr 2 \ \jf kX (n ) | 2 V (n , r 2 ) I Wk'X' {ri)\ 2 fa) 



U^'(q) = (kAk'A'|V|k'AkA') (24) 
= J dr\dr 2 y* kX (r\)y Vk (n) (25) 
x V(n,r 2 )\l4t X (r 2 )YkX'(r2) 



The matrix elements in the Left-Right (LR) basis can be 
written in the following way 



(kL\ V H \kL) = £ \L kX | 2 £ NkU^ (n k 



X'X> 



(26) 



k'X' 



where Nk is the total number of k points used in the sum. The 
density matrix operator is defined as n XX i = c^,c x where the 
label A represents the quantum numbers that labels the basis 
set. The explicit forms of the amplitude coefficients L k \ in 
the tight-binding model can be found in the appendix II. The 
diagonal n x = n xx is simply the occupation in the state A and 
implicitly implies a sum in the k points 



1 V t 

nxx> = J^h c kX c kX- 



(27) 



We get a similar expression for the other element of the 
diagonal where we only need to change the L label into R. 
The Hartree term averages to zero for the off diagonal matrix 
element. In a neutral ribbon the presence of the positive back- 
ground charge neutralizes the electrostatic Hartree-potential 



V ext (l)+V H (l) = 



(28) 



so we can focus our attention in the exchange contribution of 
the interaction. 

Now we show the matrix elements of the Fock term. For 
the diagonal and off diagonal elements we have respectively 



(kL\V F \kL) = (kL\- ^U xx '(k' -fe)(4 A ,c w )4c u ,|feL) = -£L* A L tt ,£C7i A '(fe'-fc)(n fc 



:'XX> 



(29) 



(kL\V F \kR) = {kL\-Y,Ux k \k'-k){cl k ,c V x)c\ x c kX ,\kR) = -^ (30) 
k'XX 1 XX' k' 



V 



differing only in the expansion coefficients of the two-bands The kernels in the Hartree and Fock terms of equation ( f22b . 
states. assuming a = a' and dropping the spin part, can be written as 

I 



j jXX' 



U XX \q) 



l K r o 7 1 K / \ 

-2 £ / dndr 2 \<j> (ri -Ri-T/)| v(ri - r 2 ) |0(r 2 - R, - z e ) | 2 ~ £ V eff ( Lg ) 

K i.j K i.j 

1 ^ iik'-^l^+Zj-lR.+z,,)) /, h 1 ^ i(k'-k)L«' ( 



(31) 
(32) 
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where we have defined L'j' 



Coulomb integrals V e ff 



= R; — Rj + T/ — T// and the 
can be approximated through 

the expression in equation (O 

The Hartree-Fock equations reduce to the Hubbard model 
when the interactions are reduced to onsite repulsion and the 
expressions simplify considerably. 

(kL\V H \kL) = £|L„| 2 f/((« z 



{nm}) 



(kL\ V F \kL) 



-t/£|L w | 2 («, 
I 



(33) 
(34) 



Adding both terms we obtain for the diagonal and off diagonal 
terms 

(kL\V HF \kL) = uY,\L k ,\ 2 (nia) (35) 
/ 

(L\ Vf \R) = -U"£L* kl R kl (n la ). (36) 



The last off diagonal term reduces to zero because the coeffi- 
cients L,u and have zero overlap. 



Appendix II. Analytic resolution of the tight-binding 
LR-functions 

In this appendix section we describe the exact tight binding 
wave functions for the zigzag ribbon edge states. We use peri- 
odic boundary conditions in the x-direction and closed bound- 
ary conditions in the y direction. Taking advantage of the 
translational symmetry we may Fourier transform the Hamil- 
tonian and solve a one dimensional problem with k = k x as a 
parameter. This derivation is similar to that of Malysheva et 
alM with some additional details. 

The ^-dependent one dimensional Hamiltonian of graphene 
is given by: 



■ //<■ 



( o R + Q" 

U f +e o , 



(37) 
(38) 



where Q = 2cos(/</2) and the vectors ^ have 2N coordinate, 
the N top coordinates are the A sublattice sites along the y- 
direction and the bottom N are the B sublattice sites. The 
operator R T is a N x N matrix which represents a translation 
by one unit cell in the positive y-direction and R translates in 
the opposite direction. Employing the ansatz: 



Wa 



(39) 



we may reduce Schrodinger's equation to: 
{z + Q)Wa=E\i/ b 



(40) 



In order to solve these equations with a non-trivial vector 
(Wa, Wb) we require that the matrix of coefficients have a zero 



determinant. This leads to a relation between z and the energy 
E: 



-E z + 2cos(£/2) 

l/z + 2cos(/t/2) -E 



= 



E 2 = l + (z+l/z)Q + Q 2 



(41) 



Note that for a solution z, 1/z is also a solution. We can there- 
fore, write two (unnormalized) solutions: 



Vi(z,E,k)= [ Z g Q )z" (42) 



In order to complete our description of wave functions we 
need one more condition on E or z. This condition comes 
from the edges of the ribbon. For pedagogical purposes, let 
us briefly mention the case of infinite and semi-infinite sam- 
ples. In an infinite system we require that the wave func- 
tions stay finite in the limits n — > ±°°. This leads to \z\ = 1 
which is obeyed by Bloch waves z" — e lp " with a real p being 
the momentum along the y-direction. In a semi-infinite sheet 
we require only that the wave function be finite at n — > +°° 
and a vanishing amplitude on the one edge. In order to 
satisfy the vanishing amplitude condition we need to com- 
bine the two solutions in Eq.(l42l and in order to satisfy the 
condition at infinity we require |z| < 1. Bloch waves (with 
\z\ = 1) can satisfy both conditions thus the bulk states are 
rearranged to accommodate the edge. In addition, new real 
solutions with \z\ < 1 may also satisfy the boundary condition 
on the edge with E = 0. This gives either z = — 2cos(£/2) or 
1/z = — 2cos(/<'/2) when only one of these solutions is physi- 
cal and decays away from the edge. The vanishing amplitude 
condition is trivially satisfied by setting all A (or B) site am- 
plitudes to zero. This gives the flat band states at momenta 
k which satisfy |cos(/</2)| < 1/2. This condition defines the 
region between the K and K' points. 

In the case of a ribbon with two edges we require a vanish- 
ing amplitude on both sides of the ribbon. This translates to 
^(0) = ^f(N + 1) = where N is the number of unit cells 
in the y direction (i.e., the number of chain pairs, as defined in 
the main text). Assuming that ^(n) is a linear combination 
of the solutions in Eq.(l42l we obtain the relation: 



Q 



z N -z- N 

7 N+\ _ -(N+l) ■ 



(43) 



which determines all the possible solutions of the tight binding 
model on the ribbon. Eq.(l43l has N — 2 Bloch like solutions 
and 2 real (edge like) solutions 60 . For real numbers z" produce 
exponentially decaying solutions which are localized to one of 
the edges. Similarly to the case of the semi-infinite graphene 
sheet, the real, edge-localized solutions are found in a region 
between the valleys. The term valleys usually refers to the 
nodes in the Brillouin zone where the energy vanishes. Here, 
we can only see their projection on the k axis at k = ± j^. 
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However, due to the finite size of the ribbon the edge states 
regime is smaller and extends from 2n/3a + l /W to 4n/3a — 
1 /W as will be shown shortly. 

In order to find the amplitudes of the left and right wave 
functions all we need to do is combine the two solutions in 
Eq. (l42l) in a way that the boundary conditions are satisfied. 
We may choose to satisfy the boundary condition provided by 
one of the edges while the other condition would be automati- 
cally satisfied by the correct choice of z (a solution to Eq.d43b). 
This gives: 




and the parameter a is given by the normalization: 



(44) 



a 



-2 



£[((- + QY - (z + Q)z-") 2 + E 2 (z n - z-"f] 



n=\ 



2(1-0 



2 A\-z 2N ){\+z 2 ^) 



^_ z 2(W+l))2 



(45) 



Please note that z is assumed real (but may be negative). It is 
worth mentioning that the lower energy solution (with — \E\) 
is antisymmetric about the ribbon center while the higher en- 
ergy solution is symmetric. For this (anti)symmetry to hold 
we require ¥ A (n) = ± x ¥ B (N+l - n) which is satisfied by 



To gain some intuition about interacting edge states we pre- 
fer to work with edge localized states and ^r. These are 
constructed by combining the x ¥±. Since the parameters z, Q 
and a only depend on the absolute value of E we get: 

= ^ + + ^-)^V2a((- z +Q)z"-(z + Q)z-")(^j 

= -±=(y+-y-) = V2\E\a(z n -z-")(^J (46) 

Note that the left and right wave functions are localized on 
one sublattice, determined by the boundary conditions. 

In the above discussion we left z as a parameter. However, 
the value of z is defined by the solutions to Eq|43j which un- 
fortunately does not provide a closed form. Let us define the 
inverse localization length, A such that z — exp(A). Eql43lcan 
be rewritten as 



Q 



sinh(TVA) 
sinh((/V+l)A)' 



(47) 



Then the energy of this localized solution can be simplified by 
substituting Eq.d47l) into Eq.lHTI): 



E = ± 



sinh(A) 



sinh((AT+l)A) 



(48) 
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